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COMPUTABILITY OF FILTERED QUANTITIES FOR THE 
BURGERS' EQUATION 

ERIK BURMAN* 

Abstract. We prove a priori and a posteriori error estimates for filtered solutions of the Burgers' 
equation discrctized using two different stabilized finite element metfiod, linear artificial viscosity and 

■ nonlinear artificial viscosity. We showr that the error of filtered quantities converges with an order in 
h. The constant of the estimate is independent of both the regularity of the exact solution and the 

j Reynolds number, it is moderate and primarily depending on the smoothness of the inital data. The 

convergence order and the norm in which convergence is measured are associated to the width of the 
filter. Similar estimates for the unfiltered quantities on the other hand lead to the huge exponential 
factors that are known to make estimates deteriorate for high Reynolds fiows. 

D 

■ 1. Introduction. The computability of solutions to tlie equations of fluid me- 
' chanics remains an open problem. Nonlinear effects such as turbulence and shocks 

make the stability analysis of the equations very difficult and any error estimate will 
include constants that are huge, making the estimates worthless in practice. Even in 
the simple case of the Burgers' equation it appears that the L^-error analysis for fi- 
nite element methods is incomplete, indeed the perturbation equation of the Burgers' 
equation does not seem to be stable in the L^-norm. The objective of the present 
I work is to show that for finite element discretizations of the Burgers' equation, filtered 

. quantities are provably computable, provided sufficient artificial viscosity is added to 

make the stability properties of the continuous equation carry over to the discrete 
case. The computability of the corresponding unfiltered quantities appears difficult 
to establish theoretically, due to the huge stability factors involved. 

We will mainly be interested in computations where the mesh Reynolds number 

■ is much larger than one. The Reynolds number in this context is a characteristic 
00 ■ velocity U, to be defined, times the smallest length scale of the flow (i.e. the mesh 

size h) divided by the kinematic velocity v. 

Re ^. 

Energy conservative methods, such as the standard Galerkin finite element method, 
have poor stability properties in this regime due to accumulation of fine scale en- 
ergy on the resolved coarse scales. Stabilized finite element methods, on the other 
hand, have some properties that make them particularly interesting for the solution 
of high Reynolds flow problems. It has for instance been proven that stabilized finite 

■ element methods applied to linear advection-diffusion problems have optimal conver- 
gence away from layers with a constant independent of the exact solution in the layer 
region [16l [7] . The stabilization also allows to improve the global convergence order 
for the incompressible Navier- Stokes' equations by a factor in estimates indepen- 
dent of ly, provided the solution is smooth [3]. Long time stability results have been 
obtained for the subscales in the variational multiscale method suggested in [T] . It has 
also been argued heuristically that the dissipation induced by stabilization matches 
the decay rate of the power spectrum under mesh refinement in the high Reynolds 
regime [3l[22]. Finally there are results proving that for the Burgers' equation, pro- 
vided the stabilization gives sufficient control of the discrete solution, the asymptotic 
limit of the sequence of discrete solutions is the entropy solution [T^l H] ■ 
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All these results are however of mainly qualitative nature in the nonlinear case 
since they do not give any information on what one actually can compute on a given 
scale. The error estimates are of purely academic interest since they will either ask 
for a regularity of the solution that may not hold, or they feature a multiplicative 
factor such as exp(|| Vu||looT) or exp(z^~"'^r) where usually ||Vu||l=° is very big and 
V very small in the regime that interests us. Asymptotic estimates for the hyperbolic 
Burgers equation on the other hand give no information of the accuracy of a certain 
computation with a fixed h and is therefore also of reduced interest in practice. 

We will consider two different stabilized finite clement methods in this work. First 
the classical linear artifical viscosity resulting in a first order scheme will be introduced 
to the standard Galerkin method and then a nonlinear viscosity in the spirit of that 
proposed in HI. This is a shockcapturing technique related to those proposed in [IS] 
or the entropy viscosity of |13| . 

Consider the simple model case of the Burgers' equation with periodic boundary 
conditions, on the space-time domain Q :— I x (0,T), with / := (0, 1) 

dtu + ^dxu'^ — vdxxU = in Q 

u{0,t) = for t e (0,r) (1.1) 

u{x,0) — uq{x) for X £ I. 

For the discretization we use the simplest of stabilized finite element methods, piece- 
wise affine continuous approximation with linear, or nonlinear, artificial viscosity. In 
this framework we prove estimates for a regularized error. The interest of these es- 
timates stems from the fact that the constant of the estimates are of moderate size 
and only depends on the regularity of the initial data, hence there is no dependence 
on the Reynolds number, nor of the global regularity of the solution. Error estimates 
for averaged or postprocessed quantities have already been proposed in the literature, 
see for instance [TOl [HI [H] , but to the best of our knowledge they all suffer from the 
huge exponential constants, or dependence on the Reynolds number, discussed above. 

First we discuss the L°°{0,T; L^{I)) stability of the Burgers' equation and con- 
clude that the resulting estimate includes an exponential factor of the type 

exp(||aj;u||L~T) 

reflecting a possible instability in the L^-norm. Then we introduce the finite ele- 
ment discretization and discuss the stability properties of the method in some detail. 
Finally we consider filtering of the final solution and show that the perturbation 
equation corresponding to the filtered solution has improved stability properties and 
the error may therefore be upper bounded independently of both the regularity of 
the exact solution and the physical viscosity. The estimates also give a precise rate 
of convergence in h, depending only on the filter width. Typically we will consider 
the following differential filter known from time relaxation methods (see for instance 

mm) 

-S'^Au + u = u{-,T) on I (1.2) 

with u{0) — u(l) and S the filter width. As we shall see, although || {u — Uh){-, T^)||l2(/), 
where denotes the finite element approximation of (jl.ll) . does not appear to allow 
for error estimates with moderate constants, \\{u — Uh){-,T)\\i^2(^j-^ does. Indeed, in 
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the high Reynolds number regime we prove the error estimate 



\\Sd,{u - u){T)\\l2^i) + \\{u- 7i,0(r)|lL2(,) < C{uo,T) exppoT) (^^^ ' (1.3) 

where u and Uh are the filtered exact and computational solution respectively. The 
constant in ()1.3p depend only on the intial data, the mesh geometry and the final time. 
We will use the notation Uq := sup^gj \TThUQ{x)\ and Dq := sup^gj ^dx{uo + Tr}iUo{x)) , 
where Wh denotes the L^-projection onto the finite element space. We will choose uq 
as a smooth function and by the stability of the L^-projection on regular meshes we 
have Uq < sup^g/ |wo(a;)| and Do ^ sup^.^/ \dxUo{x)\, so that estimate depending on 
Uo and Dq are indeed mesh independent. For simplicity we assume ug G C°°{I), this 
does not exclude the formation of sharp layers with gradients of order at later 
times. For fixed filter width (II. 3p results in a convergence rate of order /12. If on 
the other hand the filter width is related to the mesh size 6 ^ h" with a < i we get 
the convergence rate h ~ . The parameter S determines how strong the localization 
of the norm is. The choice 6 — 1 leads to a norm related to the iJ^^-norm and the 
choice 6 — h leads to a norm similar to the L^-norm. Clearly the estimates proposed 
here only makes sense for < a < ^. Showing that no error bounds in a norm similar 
to the case can be obtained. If the aim of the computation is to use Uh as an 
approximation of u, we observe that the L^-norm perturbation error induced by the 
filtering in the smooth portion S oi I satisfies 



Here we use the notation a <h defined by a < Cb with C a constant independent of 
h, the physical parameters and of the exact solution. 

The computational error of the filtered solution, as an approximation of u in the 
smooth part of the domain on the other hand is therefore given by 

11"" - Uh\\l^(S) < \W - ^\\l^(S) + \\u - UhWb^iS) ^5"^ + (yj^ 



and with 5 :— h°' we obtain 



Uh\\L\S) < 



It follows that the optimal convergence order away from layers is obtained for a = i, 
leading to the L^-error 0{h^) in the smooth part of the domain S. It is also worth 
noting that thanks to the maximum principle, under the reasonable assumption that 
meas(/\ S) ~ O{fj./Uo) we have the following bound for the global L^-error 

1 

\\u-Uh\\L^(i) < \\u ~ 'S/i||l2(5) + U() ^ "^^^(^^^0^) ^ <hi + {Uoma.x{^i,Uoh))i 

which implies that for Re >> 1 the global convergence order in is of order 0{h^). 
Due to the compound effect of the regularization error and the discretization error, 
most likely this estimate is sharp for ii/j, whereas Uh could provide a superior approx- 
imation of u. 

The derivation of the estimate (|1.3p uses: 
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- sharp energy stability estimates for the finite element method, 

- maximum principles for the finite element solution and its first derivative, 

- a priori stability estimates on a linearized dual problem with regularized data, 

- Galerkin orthogonality and approximability. 

The present analysis draws on previous ideas on a posteriori error analysis for hyper- 
bolic problems and computability in CFD, mainly from TT', [5^ and [Tl|. Below we 
are interested in estimates where we have precise control of all constants. This implies 
controlling the asymptotic growth of the residual and working with the exact dual 
adjoint, involving both the approximate and the exact solution. In the literature, the 
adjoint equation is often linearized around the exact solution only, eliminating the 
feedback of the stability properties of the numerical scheme in the sensitivity analysis 
and preventing the derivation of error estimates. 

The viscosity coefficient of the dual problem herein is the physical viscosity and 
not the artificial viscosity, as often has been suggested. This way the dual problem 
contains the accurate sensibility information of the physical problem. 

2. The Burgers' equation with dissipation. The wellposedness of the equa- 
tion (|l.ip for v > is well known (sec for instance [21 ) it is also known that for > 
by parabolic regularization the solution is C°°{I). This high regularity however does 
not necessarily help us when approximating the solution, since we are interested in 
computations using a mesh-size that is much larger than the viscosity and still want 
the bounds to be independent of high order Sobolev norms of the exact solutions and 
of 

2.1. L^-stabiUty of Burgers' equation. Consider a general perturbation r]{x) 
of the initial data of 

dtii + \dxU^ ~ vdxxU = in Q 

u{Q,t) = for i e (0,T) (2.1) 

Zi(x, 0) — uq{x) + rj{x) ior X £ I . 

Taking the difference of (|2.ip and (|l.ll) leads to the perturbation equation for e := u—u 
with a(u, u) := i(w + m), 

dte + dx{a{u,u)e) - vdxxG = in Q, 

e(0,t) = e(l,i) for t e (0,r) (2.2) 
e{x,Q) — r/{x) for x <E I. 
Multiplying equation (|2.2I) by e and integrating over Q leads to the energy equality 




We know that due to shock formation ~dxa{u,u) ^ v^^ . Any attempt to obtain 
control of ||e(r)||^2(j2-) in terms of the initial data will rely on Gronwall's lemma, 
leading to 

\\e{T)\\l^i,i)<CaMl^ii) 
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with the exponential factor 



Ca := exp(||aa;a(u,u)||i=e(Q)T) ~ exp(T/j/), 

which is exactly the bad factor discussed in the introduction. It tells us that we have 
stability (and hence computability) only up to the formation of shocks. Using this 
type of argument in the analysis of the finite element method leads to error estimates 
of the type derived in [6], useful only for solutions with moderate gradients. 

2.2. MELximum principles for Burgers' equation. It is well known that the 
equation (jl.ip satisfies a maximum principle on the form: 

sup t)| < sup |uo(a;)|. (2-3) 

{xJ)GQ xei 

This follows using the smoothness of the solution u and the standard technique of 
assuming a local extremum in a point and showing that dtu has the appropriate sign 
in that point (or in the hyperbolic case, using the method of characteristics). For our 
purposes we also need some precise information on the derivative. Since the solution 
of (|1.1[) is smooth we may derive the equation in space to obtain the following equation 
for the space derivative w := dxu: 

dtw + udxW — vdxxW = —vo^ in Q 

w(0,t) = w(l,t) for t € (0,T) (2.4) 

^(x, 0) — dxUt){x) for X <E I. 

Assuming that w takes its maximum in some point x (z I and noting that dxw{x) — 
and dxxw{x) < it follows that dtw < Q aX the maximum and we deduce the bound: 

max dxU < uY&yidxUQ. (2-5) 

It follows by the smoothness of the initial data that the space derivative is bounded 
above for all times. 

3. Artificial viscosity finite element method. Discretize the interval / with 
N elements and let the local mesh-size be defined by h := 1/N. We denote the 
computational nodes by := i h, i = 0, . . . , N , defining the elements Ij := [xj,Xj+i]^ 
j = 0, . . . , N — 1, and the standard nodal basis functions {vi}^Q, such that Vi{xj) = 
5ij , with 6ij the Kronecker delta. To impose periodic boundary conditions we identify 
the node xq with xn and define the corresponding basis function vqn '■ ixo,xi) U 
{xn-i,xn) I— K by uq on {xo,xi) and by vn on {xn-i,xn)- This basis function then 
replaces vo and vn, leading to a total of N degrees of freedom. For simplicity we use 
the notation vq for the basis function vqx- The finite element space is given by 

{N-l 
UiVi, where {ui}^o^ ^ * 
4=0 

We define the standard inner product on X C / by 



{vh,Wh)x ■ 




VhWh Ax. 



The discrete form corresponding to mass-lumping reads 

N-l 

{vh,Wh)h ■= ^ Vh{xi)wh{xi)h. 

i=0 

The associated norms are defined by := (w, for all w e L'^{X), if X coincides 

with / the subscript is dropped, and H'yhH/i ■— {vh,Vh)f^ for all Vh G Vh- Note that, 
by norm equivalence on discrete spaces, for all Vh G Vh there holds 

\\Vh\\h < \\vh\\ < \\vh\\h- 

Using the above notation the artificial viscosity finite element space semi-discretization 
of (jl.ip writes, given uq g C°°(ri) find Uh{t) G Vh such that {uh{0),Vh)i = {uo,Vh)i 
and 

{dtUh,Vh)h + (^dx^,Vh^ + {C'dxUh,dj;Vh)i = 0, for all Vh € Vh and t > 0, (3.1) 

where we propose two different forms of 

1. linear artificial viscosity: 

i) := max(C/o/i/2,zy); (3.2) 

2. nonlinear artificial viscosity: 
Let < e and 

/ M 1|| II \[dxUh]\x\ 

Muh)\h — o "ft L-(70 max — — , 3.3 

2 ' xe{xi,xi+i} 2{\dxUh\}\x + e 

where [9a:M;i]|£Ci denotes the jump of dxUh over the node Xi and {|9a;U?i|}Ui 
denotes the average of IcJicW/il over Xi. If e = and {|9a;Uh|}|a;i = we replace 
the quotient |[9xM;i]UJ/{|9a;M/i|}Ui by zero. 
Further let 

1 if dxUhlu > 0, dxUhlu > dxUhlu^^ > 
^{uh)\h ■= { and dxUh\i^ > dxUh\i^_i > 
otherwise 



Muh)\h c("/i)l/. max! i'o|/,_i-^ . — ,J^o|/i+i-^ . — • (3.4) 

Finally define: 

v{uh)\i, := max(i^, /i(i'o|/i + J^i|/J). (3.5) 

The rationale for the nonlinear viscosity is to add first order viscosity at local extrema 
of the solution Uh so that (|2.3p holds also for the discrete solution and enough viscosity 
at positive extrema of dxUh, making ()2.5p carry over to the discrete setting. The most 
important term is I'o, ensuring the discrete maximum principle. The other part vi is 
merely a correction that ensures that the viscosity at local maxima of dxUh dominates 
that of the surrounding elements. Indeed the role of the function ^{uh) is to act as 
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an indicator function for the elements where the local maxima of d^Uh are taken and 
modify the viscosity there. By construction, if ^(w/i) = 1 in one element it must be 
zero in the neighbouring elements. Note that formally i'{uh) « 0(max(L/o^'^^^)j 
in the smooth part of the solution, so that in principle we can expect higher order 
convergence away from local extrema. Note that the perturbation vi close to extrema 
of the gradient of the solution has no impact on the formal order. The high order 
convergence properties and the effect of the regularization parameter e will be explored 
in the numerical section. 

Remark 1. Note that in the linear case the viscosity v may he written as 

V :— i/max(l, i?e;i), with Rch ■= — — 

21/ 

reflecting that in the high Reynolds number regime the viscosity is increased artificially 
to be order h. 

3.1. Existence of solution to the semidiscretized system. For the linear 
method existence and uniqueness of solutions of (|3.1|) follows using standard methods. 
First we observe that the nonlinear function of the dynamical system is locally Lip- 
schitz and then anticipating the global upper bound p.6p we conclude that a global 
solution exists and is unique. 

The nonlinear method obtained when (|3.1|) is used with the viscosity p.5|) results 
in a dynamical system with discontinuous righthand side (even for e > the contri- 
bution from vi introduces discontinuities). Existence of solutions to (13. ip with the 
nonlinear viscosity p.Sp is obtained using Filippov theory [11] . Anticipating the re- 
sults of the next section, we may conclude that a solution exists, since by the discrete 
maximum principle p.6p . for fixed h, there exists M > such that 



< M 



and hence \dtUi^{x,t)\ < M for all {x,t) G Q. The question of uniqueness is more 
involved, but using the following argument we conjecture that the solution is forward 
unique. Observe that in a given pair of neighbouring cells A/ :— U li, the 
system can either have ^ = in both cells, or have ^ = in one of the cells and 
have ^ = 1 in the other. It follows that there are three states possible in A/ across 
the boundaries of which the nonlinear function has a jump. Denote the interiors of 
these configurations by 5*1 {^{uh)\Ai = 0), 5*2 = 0, = 1) and 5*3 

— 1; — 0) the boundary between Si and Sj may be denoted 

Yiij. By construction (see the bound p.Zp ') the nonlinearity will drive any solution 
in S2 or S3 towards 5*1 and no solution whose restriction to A/ is in Si can enter 
S2, or S3. This rules out the so called repulsive sliding mode that is known to cause 
nonuniqueness of the solution. We will not explore these issues further here but refer 
the interested reader to [Sj. Typically in practice the system p.ip will be discretized 
in time using an explicit time stepping scheme which, by definition, will produce a 
unique discrete solution. In the following we will prove that any solution to (j3.ip will 
satisfy certain uniform bounds and converge to the exact solution at a certain rate. 

3.2. MELximum principles for the discrete solution. Maximum principles 
give local estimates of the behavior of the solution and they rarely carry over to the 
discrete method. There are however stabilized methods that are specially designed to 
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make a discrete maximum principle hold, see for instance |23[ [5] for linear convection- 
diffusion problems and |19[ |1] for maximum principle satisfying finite element methods 
for conservation laws. 

In [4] it was shown that the discrete equivalent of (|2.3[) . together with energy 
stability of the discrete solution is sufficient to prove the convergence of the approx- 
imation sequence to the entropy solution. For our purposes herein however it is not 
sufficient, but we also need to prove a discrete equivalent of the bound (|2.5|) on the 
gradient. We collect the monotonicity results we need in the following lemma. For 
clarity of the exposition we first give the proofs in the case e and then discuss 
how the regularization modifies the bounds. 

Lemma 3.1. Let he the solution of (|3.1[) either using the linear viscosity (|3.2p 
or the nonlinear viscosity p.5p and e = 0. Then the following hounds hold: 

sup \uh{x,t)\ <Uq <m&yi\uQ{x)\, (3.6) 
(a:,*)eQ ""^^ 

sup dxUh{x,t) <Tn&yidxUh{x,Q) <Tnsx\dxUQ{x)\. (3-7) 
(x,t)GQ ""^^ ^s/ 

Proof. The proof of (I3.6P is an immediate consequence of the fact that the space 
discretization has the DMP-property introduced in [5 . 

For the case of linear artificial viscosity, first assume that for some time t* there 
holds maxa;g/ < U^, then show that this implies 

maxmax |u/i(a;, < J7o (3-8) 

t>t* xG/ 

and conclude noting that the assumed inequality holds for = 0, since 

max|u,i(x,0)| =: Uq. 

xG/ 

First we compute 



f h \ 

/ UhdxUhVi dx = -dxUh\]., + -^dxUh]]. 
J I 3*^6' 



+ ^Uh{x.i^i)dxUh\u_^ + ^Uh{xi)dxUh\u (3.9) 
Assuming that Uh has a local max in Xi aX t ~ t* it follows that, for linear viscosity 



hdtUh{xi,t*) = -{uhdxUh,Vi)j - [OdxUh,dxVi)i 

< (Ua - max(|u,i(a;i_i)|, \uhixi)\)) {\dxUh\i._-^ \ + |9^u/i|/J) < 0. 

It follows that dtUfi{xi, t*) < and hence the local maximum can not grow. The case 
of a local minimum is similar. 

For nonlinear viscosity on the form (j3.3l) . since at a local maximum = 
2{\dxUh}, we deduce that 

^'("/i)|/i = \\uh\\L-^(ii)h/2 = ma.x{\u{xi)\, \u{xi+i)\)h/2 

and the same conclusion follows. 
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We will show (j3.7p by first proving that the maximum gradient must be decreasing, 
and then applying the stability of the projection. Hence it is sufficient to prove 
that dt maxi dxUh\ii < in Q to conclude. 

We first give the proof for the linear artificial viscosity. We will prove that the 
discrete gradient is bounded by the gradient of the discrete initial data. 

sup dxUh{x,t) <s\y(>dxTThUa- 

Starting from p.ip . we let li be any element where dxUh\ii has a local maximum, in 
the sense dxUh\u > dxUh\u±i > 0. 

I j■x^+2 

dtdxUh\u = -T UhdxUh{vi+i - Vi) dx- 

vdxUhdxivi+i - Vi) dx ^ Ti + T2. 



h 

Decomposing the integral Ti on the contributions from Vi and w^+i we have using 
p.9p and some minor manipulations 

\h{dxUh\u_,)^ ^'^KdxUhU,)^ 

]^Uh{xi){dxUh\i^ - dxUh\u_^) - ^ 



Ti = --h{dxUh\u_,f - ^KdxUh\i:)' ~ -h(dxUh\u^,f 

^Uh{xi){dxUh\u - dxUh\u_^) - ^Uh{xi+i){dxUh\ii^i - dxUhlu). 



Since the derivative takes its max value in li we have for 
I rXi+2 



h 



udxUhdxivi+i - Vi) dx 
= h^^ {{v{uh)dxUh)\u_^ - 2{i'{uh)dxUh)\u + {i'{uh)dxUh\u^i)) < 0. 
Collecting the above expressions and using that i> > ^Uoh, and we obtain 

Ti+T2< -^KdxUhlu^,)^ - ^HdxUhlu)^ - ^h{dxUh\u_J^ 
- (C^o + Uh{xi)) {dxUhlii - dxUh\ii_J 

^ V V ' 

>0 >0 

- (C^o - Uh{xi+i)) {dxUh]ii - cJ^jW/il/.+i) < 0. 

^ V V ' 

>0 >0 

This proves that dtTaaxi OxUhli^ < and therefore the maximum space derivative is 
alway decreasing. The global upper bound (j3.7l) is immediate by the stability of the 
i^-projection. Note that by the non-monotonicity of the L^-projection there is now 
an absolute value on the derivative of the initial data. 

In the case of the nonlinear viscosity given by p.Sp we first show that the time 
derivative of the gradient must be negative in cells with £,{uh) = 1. Then we show that 
any adjacent element, cannot grow either due to the design of the nonlinear switch. 
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In this case have after integration 

Ti+T2 = -lh{d.,uu\u^,f - lh{d^UH\u? - lKd,Uh\u_,f 
6 3 6 

- h^^{h'h{uh)\iidxUh\h - Vh{uh)\u_^dxUh\u_^) 

- h^^{vh{uh)\iidxUh\u - C'h{uh)\u+idxUh\u+^)) 

- ]^Uh{xi){dxUh\u - dxUh\u_^) + ]^Uh{xi+i)[dxUh\u - dxUh\u^^). (3.10) 

First observe that the case where either Xi or x^+i is a local extremum can be excluded, 
since then vi^\i- — ^-nd by observing the sign of the contribution of the 

derivative from the neighbouring cells in the viscosity terms of line two and three of 
p.lOp . The other terms are controlled as in the linear theory with minor modifications. 
Since the gradient has a local max in 7^, in the sense that £,{uh)\ii = 1 there holds 

viuh)\ii = (max(i'o("h)|/i_i^:^ i — ,va\Uh)\h+t^ j — ) + vo{Uh)\h) 



dxUh\h ' dxUh 



li 



and since by construction C(''^/i)l/i±i = 0; have in the neighbouring cells, 

^'("/i)|/.±i = ^^o('i/i)|/i±i- 

Using the values of i) the contribution from the viscous part of the differential operator 
may be bounded as 

- h~'^{vh{uh)\hdxUh\ii - i'h{uh)\h-xdxUh\ii_i) 

- h^^{vh{uh)\hdxUh\ii - yh{uh)\h+tdxUh\h+i)) 

< -2vn{uh)\i.dxUh\u. (3.11) 

For the two last terms in the right hand side of p.lOp we note that, assuming first 
Uh{.Xi+i) > Uft(xi) > 0, 

- T:Uh{xi){dxUh\i^ - dxUh\ii_:^) +l:Uh{xi+i){dxUh\i^ - dxUh\i,+i) 
Z V ^ / z 

>o 

< Muh)\hi\dxUh\u\ + \dxUh\u^A) (3.12) 

Collecting p. lip and p.l2p . p.lOp can be upper bounded in the following fashion, 
recalling that d^Uhlu > d^Uhlu^-, > 0: 

Ti+T2< -lh(dxUh\u^,f - h{dxUh\uf - h{dxUh\u_,f 
b 6 b 

-2i^o{uh)\hdxUh\i^+i^o{uh)\u{\dxUh\ii\ + \dxUh\u^i\) < 0. (3.13) 

The case > Uh{xi+i) > Uh{xi) is similar observing that in that case the last term in 
the right hand side of p.lOp is negative. In case Uh{xi) < < Uh{xi+i) we observe 
that only the treatment of the first contribution of the last line of p.lOp must be 
modified. We note that 

-]^Uh{xi){dxUh\i^ - dxUh\i,_i) < ^h{dxUh\u)^ 
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and that this term is cancehed by the second term of the right hand side in the first 
hne of (IXTU)) . 

FinaUy we must check that the gradient can not increase in any portion of the 
domain where the gradient is constant at the maximum value over several elements 
from Im to /„ at some time t* . First note that for all elements Im+2, In-2 the 
derivative is decreasing, since only the first three terms of the right hand side of (|3.10p 
are non-zero. By construction C(^^?i)|/„ = 1 E^nd hence the derivative is decreasing in 
/„. As a consequence the derivative in is either decreasing at the time t* or will 
have = 1 at t* + £ for all £ > and hence be non-increasing. Similarly for 

Im and Im+ii the derivatives can not grow at the same rate in both cells since then 
Cl/m+i = ^ a,t t* + e for all £ > 0, since its right hand side neighbour has decreasing 
space derivative. If, on the other hand the time derivative of the derivative is the 
largest in at time t = t* then ^|/^ = 1 at i* -I- £ for all £ > 0, and hence the 
derivative can not grow. This case of several adjacent cells over which the gradient is 
constant is what can give rise to the so called attractive sliding mode in the Filippov 
theory. □ 

Remark 2. First observe that if vi — it is not difficult to find Uh for which 
dtTaaxi dxUh\ii > 0, so the above technique of proof requires the contribution from vi, 
whether or not it is really necessary in practice remains unclear. 

The form of vi can be simplified. Lemma \3.1\ also holds for 

viiuh)]!, inax(^t'o|7,_i,t'o|/i+i), (3.14) 

Since —fT-, — < 1- 

It was shown in [4J that a consequence of the bound p.ip is that the total variation 
of u/i diminuishes. We recall the result without proof. 

Corollary 3.2. Let be the solution of p.ip . then there holds, for all t >0 

TV{uh{-,t)) -.^ j\dxUh{-,t)\ dx<TV{uh{-M- 

3.2.1. The effect of non-zero regularization parameter e. In practice it 
may be practical to use a value on e that is related to the mesh size, in particular if 
implicit solvers are used. This will result in a modification of the upper bounds p.6p 
and (13. 7|) . but as we show below, the maximum principles can only be violated by an 
0{e). 

Proposition 3.3. Let < eT < I in p.3p . let Uh be the solution of p.ip . then 
there holds 

lk/.(-,r)|U^(Q) < (l + 6T)f/o (3.15) 

and 

max dxUh(x,t) < i<[ia:x.dxUh(x,0) + Uo{l + eT)eT. (3.16) 

Proof. Assume for simplicity that u/i > and that u/i takes a global (positive) 
maximum in Xi that will grow with the maximum rate throughout the computation. 
Introduce the notation 



gi ■■= \dxUh\u_i\, g2 ■= \dxUh\u\- 
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Recall that at a local maximum of M/j, ^(m/j) — and therefore ^'(u/i) = majc{i>,i>o{uh))- 
Assume that the maximum is taken for VQ{uh). Then by p.9p 

dtUh{xi) < ^Uh{xi){gi +52) - h^^i'h{uh)\h-igi - h~^Vh{uh)\hg2) 

^1 . w , (.91 +92? .1 , . 

2 31 + 52 + e 2 

By Gronwall's lemma it follows that 

Since e'^ < 1 + x/{l — x) we conclude, using the assumption that eT < 1, 

Uh{x^,T) < Uo{l + eT). 

To obtain the inequality p.l6p we reason in a similar fashion starting from the equa- 
tion p.lOp and using p.lSp . The regularization only comes into effect at the step 
p.l2p and we observe that in this case 



^Uh{xi){d^Uh\u - dxUh\u_^) + ]^Uh{xi+i){d,^Uh\u - dxUh\u+^) 

< i^o{uh)\ii{\dxUh\\ii + \dxUh\\u+i +e). (3.17) 



This then leads to the bound 



1 2 1 1 

dtdxUhlu < ~Qh{dxUh\u_^i)'^ - -^h{dxUh\ii)'^ - -^h{dxUh\i,_J'^ + -Uh{xi+i)€. 



Integrating in time shows that 

max dxUh < max9a;M/i(a;, 0) + Uo{l + eT)eT. 

□ 

These perturbations of the discrete maximum principles then modifies the result p.2p 
in a nontrivial way. This will not be explored herein and we will only consider the 
case e = for the stability estimates below. 

3.3. Energy stability. Our estimates rely on stability of the numerical scheme 
and regularity of the dual perturbation equation. We need to control certain Sobolev 
norms of the discrete solution in energy type estimates similar to that of the continuous 
problem. The proof of the below estimates can be simplified in the linear case and 
the inverse estimate on the i°°-norm that is only valid in one dimension can then be 
avoided. Here we only give the proof valid both in the linear and in the nonlinear 
case. 

Lemma 3.4. The solution of the formulation p.l[) with either the linear 
artificial viscosity given by (13. 2p or the nonlinear one of p.Sp with e = 0, satisfies 
the upper bounds 

\\uh{T)\\ + \\0idxUh\\Q < \\uo\\ (3.18) 
Wdtuhh < ({/oTH-3 + i.-2)\\dxUo\\. (3.19) 
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Proof. The estimate p. 181) is immediate by taking Vh — uj^ and noticing, by inte- 
gration by parts and the periodic boundary conditions, that the nonhnear transport 
term vanishes. By norm equivalence and the stabihty of the i^-projection 

W{T)\\ < \\uh{T)\\h and ||u/,(0)||h < ||mo||. 
The second estimate fohows by taking Vh = dtu^ to obtain 

\\dtUh\\l dt ^ - {uhdxUh,dtUh)i dt- I {vdxUh,dxdtUh)i dt. (3.20) 

First note that by Coroharv 13.21 and p.6p we have, since TV(uh) < nieas(/)^ ||(9^w/j|| , 

{uhd.^uh,dtuh)i dt < UoTViuhi;0)) [ \\dtUhi;t)\\L^(^i) dt 

Jo 

< Uo\\d^Uh{-,0)\\Thri\\dtUh\\Q (3.21) 
For the last term in the right hand side of p.20p we observe that, 

i-T 

{vdxUh,d^dtUh)i dt 

T pT 

{ma.x{0,i) - i^)dxUh,dxdtUh)i dt + / {vdxUh,dxdtUh)i dt 

Jo 

<TV{uh{-M pd:,dtUh\\L^ii) dt+'^\\vd^Uh{-,T)f-]^\\vd^Uh{-M\'' 

< Un\\d^Uh{-,0)\\h~^Ti\\dtUh\\LHQ) 

+ ^\\,yidxUhi;T)\\^ - ^\\iyid,Uh{;OW. (3.22) 

Hence by applying (|3.21l) and p.22p in the right hand side of p.20p and norm equiv- 
alence in the left hand side of p. 201) we obtain the bound 



119*^^11$ < c, I WdtUhWi dt 



The conclusion is immediate. □ 

4. The linearized dual adjoint. We introduce the linearized adjoint problem 

-9tV9 -I- a(u, uy^d-j,ip - vdxx^p = in Q, 

ip{0,t) = (^(l,t) fori G (0,r], (4.1) 

ip{x,T) — ij]{x) for a; e /, 

where a(u,Uh) '■— {u + Uh)/2. The rationale for the dual adjoint is the following 
derivation of a perturbation equation for the functional of the error |(e(T), where 

13 



e(T) := u{T) - Uh{T). 



|(e(r),V')/| = |(e(r),V')/ + / {e,-dtif + a{u,Uh)da:ip-vda:xV)idt\ 

Jo 

= |(e(0),^(0))/+ / idte + dx{a{u,Uh)e),ip)i dt+ [ {i^d^e, dxip)i dt\ 
Jo Jo 

= \{e{0),ip{0))i ~ {dtUh + UhdxUh,ip)i dt- {vdxUh,dxip)i dt\. (4.2) 



This relation connects the error to the computational residual weighted with the so- 
lution to the adjoint problem and can lead both to a posteriori error estimates and to 
a priori error estimates, provided we have sufficient information on the stability prop- 
erties of the numerical discretization methods and of the dual problem. The a poste- 
riori error estimate uses techniques similar to the now classical dual weighted residual 
method, however in our case we can estimate the dual weights analytically, accounting 
for perturbations, both in the discrete and the continuous solution. Combining the a 
posteriori bounds with strong stability properties of the numerical method, leads to a 
priori upper bounds of the a posteriori quantities, showing that these must converge 
and in consequence that the error goes to zero. Before proceeding with this analysis 
we derive an a priori estimate for the derivatives of the dual adjoint (|4.1|) . 



4.1. Wellposedness and stability. Since a{u,Uh) & W^'°°{I) the problem 
(|4.ip has a unique solution and one may show that it satisfies the maximum principle 

max |(/3(a;,i)| < max I ■(/'(a;) I . 

(x,t)GQ x£l 

The following stability estimate follows easily by standard energy methods 
Lemma 4.1. Let tp be the solution to (14. ip then there holds 

sup \\dM-,t)f + '^\\dxxip\\l<exp{DoT)\\dx^f. (4.3) 
te(o,T) 

Proof. Multiply the equation (14.11) by —dxxf and integrate over (t,T) 
\\dxV{-,t)f + 2 u\\dxx^\\^ dt ^ \\dx^^ - 2 {a{u,Uh)dx^,dxxV) dt. (4.4) 
Note that, by an integration by parts and the maximum principles (|2.3p and p.7p 



(a(u,Uh)dx^,dxxV)i dt = - / {dxa{u,Uh)dx^p,dx^) dt 



DoJ^ Wdx^W^ dt. 



> - 



The result in L°°{0, T; L'^{I)) follows from the Gronwall's lemma and taking the supre- 
mum over t e (Oj^) of the resulting expression 

\\dxV{-,t)\\^ <exp{Dot)\\dxi^\\^. 
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The result for the second derivatives then follows by using this expression to bound 
the right hand side of (|4.4I) 



2/ '^\\d.j;xf\\'^dt<\\dxij\\'^+ {dxa{u,Uh)dx(p,dxf) dt 
Jq Jo 



< (l + DoJ expiDot) di) = exppo^)!! 9,^11'. 

□ 

5. Error estimates for filtered quantities. We will consider the differential 
filter defined in ()1.2p . where S denotes a filter width to be specified. The norm 
associated to the differential filter is given by 

\\\u\\\s:^i\\Sd,u\\' + \\uf)^. 

We introduce the filtered error e := u — Uh, where u and Uh denote the filtered exact 
and approximate solutions respectively obtained by solving (|1.2p with u and as 
right hand side. The analysis uses the stability properties of the adjoint perturbation 
equation (Lemma 14. ip and the stability properties of the discrete problem (Lemma 
13. 4[) to derive first a posteriori error bounds for the filtered quantities and then a 
priori bounds by upper bounding the a posteriori residuals, by a priori quantities. 

Theorem 5.1. Let u be the solution of (jl.ip . be the solution of (j3.ip . Then 
the following holds: 

• A posteriori upper bound 

\\\e{T)\\\s < eMDoT) ( ^) (hi\\{u - u,,)iO)\\ + inf \\vh + Uhd,Uh\\ dt 

\o'^ J ^ Jo vheVh 

+ [ Wd^dtUhW dt+ [ 
Jo Jo 



i-T rT ^ 

I max(0, — i^) 2 dt 



fT 1 

+h(^j^ i^muhWii dtyy (5.1) 



where 



\\[dxUh]\\N ■= '^{dxUh{xi)\i^^i - dj;Uh{xi)\i^ 



i=0 



with In identified ith Iq by periodicity. 
A priori upper bound 



\e\\\s < 



cxp{DoT) (-^y (^(^hi+ujVT')\\uo\\ + {TUo + hK^)\\dxU„\\ 



(5.2) 

Proof. Let ip = e{T) in the definition (|4.ip of the dual adjoint problem. Using 
the design of the dual problem and the definition of e we have by the relation (|4.2p 

|||e(T)|||^ = {6dxe{T),dxe{T))i + (e(r), g(r)), = (e(r), g(T)), 

= {e{0),ip{0))i - {dtUh + UhdxUh,(p)i dt- {vdxUh,dx^)i dt. 
Jo Jo 
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Taking — i^hVi with Tt^ denoting the standard L^-projection, in p.ip and adding 
to the above expression yields 



Pllla = (e(0),',i'(0))/ - / {dtUh +UhdxUh,(p)i dt - {vd,^Uh.,dx'-p)i dt 

Jo JQ 

Ja Jo 

= (e(0),<^(0))/- / {dtUh + UhdxUh,^-T^h^)i dt 



To 



Ti 

{{dtUh,iTh(ph)i - idtUh,nh(ph)h) dt- / {i^dxUh,dx{ip - nhip))j dt 



T2 T3 

+ / (niax(0, v{uh) - v)dxUh, dxT^h^)i dt . 



Now consider the terms Tq to r4 term by term. First use the orthogonahty (cq, Vh)i 
for all Vh &Vh, 

To = (e(0),i^- 7r,i(p)/ < ||eo||/i sup \\dxip{t)\\. 

te{o,T) 

Similarly for all Wh G Vh there holds 



T 

{wh + UhdxUh, ifi - ■nh'p)i dt 



rT 

< I \\wh + Uhd^UhW dt sup \\{(p - Trhip){t)\\ 

t6(0,T) 



and hence 



Ti<h inf \\vh + UhdxUhW dt sup 

Jo ^'h<^Vh t£(0,T) 

Let Ifi denote the standard Lagrange interpolant. By the definition of the discrete 
i^-inner product (•, ■)h there holds 



T2 = - y J (dtUhiTh^ - Ih{dtUhT^h^)) dxdt 

< / h'^ldxdtUhdxTTh'fl dxdt < / h'^\\dxdtUh\\\\dxTTh'p\\ dt 
Jo J I Jo 



< 



T 

2 



h \\dxdtUh\\ dt sup \\dxipit)\\ 
te(o,T) 
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For we have after an integration by parts and using a trace inequality followed by 
approximation 

T N-l 



< 



T 

i^\\[dxUh]\\N{h~'^\\f - nhfW + h2\\da;{(p - TTh^)\\) dt 







3 ,, 1 



Finally the non-consistent artificial viscosity term is controlled using the Cauchy- 
Schwarz inequality and the i?^-stability of the i^-projcction |19a;7r/i(^|| < H^xV?] 

rT 

Ti < ma.xi>\i. \\ ma,x{0, i) - iy)idxUh\\ dt sup 
' Jo te(o,T) 

<{Uah)^ / \\ rRax{0, D ~ iy)^d^Uh\\ dt sup \\d^f{t)\\. 
JQ te(o,T) 

Collecting the bounds for Tq — T4 and using the stability estimate (|4.3p we have 



ll|e(T)|||^<exp(i^or) (^^ j ' (/i^ ||e(0)|| + /i^ Jn^^Jv^, + uAu^W dt 

+ C/g^ / \\uisyi{Q.v - v)^dxUh\\ dt + h"^ / \\dxdtUh\\ dt 
Jo Jo 

1 

T 



+ hi i.\\[dxu^]\\l dt\ )mt)\\\s (5.3) 





from which (|5.ip follows. 

The a priori error estimate now follows by using discrete stability to bound the 
residuals. Since we do not assume any regularity of the exact solution we can not as- 
sume that any stronger bounds hold. Note that by a well known discrete interpolation 
estimate there holds, 

inf \\h^Vh - UhdxUh)\\ < h\\[uhdxUh\\\ n ■ 

Vh&Vh 

Using that 

h[uhdxUh]\x, < h\uhixi)\ (^ 2{\d^uh\} + e ) '"^^ (l^^^'^IUx-i + \dxUh\\u + e) 

< {i'iuh)dxUh)\u_i + {i'{uh)dxUh)\u + hUoe 

we may deduce 

inf \\h^{vh -UhdxUh)\\ < U^WO^dxUhW +h^Uoe. 
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In the linear case the inequality is trivial by taking = and using (13.61) . 

Then use the Cauchy-Schwarz inequality in time for the two terms of the second 
line of (j5.3l) . an inverse inequality for the second term in the second line and a trace 
inequality for the last term of (|5.3p . 

ll|e(T)|||, < eMDoT) (^^) ' (/i^||e(0)|| + {h^ + VT)\\i.i d,u,,\\Q 

+ h^VT\\dtUh\\Q + T^hiUoey 
We conclude by applying the stability estimates and p.lSp and p.l9p leading to 

lllelll, < eMDoT) (^^^ ' [{hi + U,^Vt)\\uo\\ 

+ iTUo + hh-2)\\d^uo\\+T-2h-2Uoey 

□ 

Observe that the above estimate is independent both of the regularity of the exact 
solution and of the flow regime. 

Remark 3. Note that if instead the initial data u/i(x,0) is chosen as the nodal in- 
terpolant of uq, ThUQ, we may define Uq — max^^g/ |uo(2;)| and Dq = maxxei dxUo{x). 
On the other hand we can no longer use -orthogonality in the upper bound for Tq. 
It appears that in that case we must use the maximum principle of the dual problem 
to obtain 

To = (e(0),^)/ < ||e(0)|Ui(,)||^(0)|U=.(,) 

< l|e(0)|Ui(/)||vI'|U~(/) <ri||e(0)|Ui(,)|p|||^ 

^ ||a>o||LM/)IPIll5, s = 1,2. 

The global convergence order will be the same, but it appears that the error contribution 
from the initial data will be larger and the factor ||9|mo||li(/) must be added to the 
right hand side of (15.21) . Another downside to this approach is that it only works in 
one space dimension, whereas before only the energy stability estimates of Lemma \3.4\ 
used one dimensional inverse inequalities. 

Remark 4- We have kept the dependence on the regularization parameter e in the 
above proof. This shows the effect of regularization on the computational error under 
the assumption that p.2p still holds under regularization. 

Remark 5. Since all these estimates are independent of v they are also valid for 
the purely hyperbolic case, with u the entropy solution. 

6. Numerical examples. In this section we will study two numerical examples 
computed with v = Q. We first consider a problem with smooth initial data 

""0 = ^(cos(7ra;) + 1). 

We compute the solution at T = 0.5, before shock formation and compute the exact 
solution on a mesh with 6400 mesh points using fixed point iteration. The intial 
data and the final solutions are given in Figure 16.11 In Table 16.11 errors in several 
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Fig. 6.1. Left smooth initial condition; right solution at T = 0.5 



N 


Wu-uhWmi) 




lllSllli 


IIISIIU 


100 


2.5- IQ-^ 


3.6 • 10-^ 


3.0- 10-* 


3.2 • 10-=* 


200 


6.7- 10--* (1.9) 


1.0 • 10-^ (1.8) 


7.0-10-^ (2.1) 


9.5 • 10-* (1.8) 


400 


1.8 • 10-* (1.9) 


3.0 • IQ-* (1.7) 


1.7-10-^ (2.0) 


2.9 • 10-^ (1.7) 


800 


4.6 • 10-^ (2.0) 


8.9 • 10~^ (1.8) 


4.2- 10-« (2.0) 


8.7-10-'^ (1.7) 



Table 6.1 
€ = 0(10""'^®), smooth solution 



different norms is presented on four consequtive meshes. Here e = to machine 
precision. Experimental convergence rates are given in parenthesis. In the foUowing 
table (Table we present the same results for e = h. The results are similar, with 
the difference that the regularized method gives second order convergence in all norms, 
whereas the one without regularistion exhibits a slight reduction in the order in the 
L^-norm. This is not surprising since formally the order of the method is and the 
unregularized method adds nonconsistent first order viscosity at local extrema. 

Now we consider a problem with non-smooth solution. The initial data and final 
time exact solution is given in Figure 16.21 We compute the solution at T = 0.5 
when the shock has formed. The exact solution is computed using the method of 
characteristics on a mesh with 12800 elements. We present tables with the same 
errors as in the previous case for the method without (Table 16. 3|) and with (Table 
16. 4p regularization. For this case there is even less difference between the two cases. 
We observe first order convergence for the L^-error and the iJ-^-norm error and 1/2- 
order convergence in the L^-norm. The computations clearly show how the weaker 
norm behaves either as an L^-norm for (5 = 1 or an L^-norm for S = h. Intermediate 
values of d appears to interpolate between these two norms. This indicates that the 
principle of our estimate, with the order depending on how S is chosen with respect 
to h is correct. A superconvergence of approximately half an order is observed for 
all the computations with the nonsmooth solution, compared to what is predicted by 
theory. 

7. Discussion. Above we proved that the error in the filtered finite element 
solution to the Burgers' equation allows for convergence proofs with order 0{h'^), e G 
(0, 5], for filter width S = and moderate constants, even though the analysis 

for the unfiltered quantity leads to error bounds with huge exponential constants. 
The results hold both for linear artificial viscosity and nonlinear viscosity of shock 
capturing type. To our best knowledge this is the first time it has been rigorously 
proven that stabilized finite element methods combined with filtering of the final result 
improves computability. 
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N 


Wu-uhWlhi) 


\\U - M/i|U2(j) 






100 


1.9- 10"^ 


3.0-10-^ 


2.3 - 10-4 


2.6-10-^ 


200 


4.7-10-'' (2.0) 


7.7-10-" (2.0) 


5.5- 10-" (2.1) 


7.1 - 10-4 (1.9) 


400 


1.2-10-4 (2.0) 


2.1 - 10-4 (1.9) 


1.3- 10-" (2.1) 


1.9-10-4 (1.8) 


800 


3.0 - 10-" (2.0) 


5.5-10-" (1.9) 


3.3- 10-« (2.0) 


5.3- 10-" (1.8) 



Table 6.2 
e = h, smooth solution 



Fig. 6.2. Left nonsmooth initial condition; right solution at T = 0.5 
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N 


\\u-Uh\\L^(I) 


\\u-Uh\\L^I) 






100 


0.036 


0.071 


6.4 • 10-'^ 


0.038 


200 


0.018 (1.0) 


0.049 (0.5) 


3.2-10"^ (1.0) 


0.024 (0.7) 


400 


9.4 • 10-^ (0.9) 


0.034 (0.5) 


1.6-10-^ (1.0) 


0.016 (0.6) 


800 


4.7-10-^ (1.0) 


0.023 (0.6) 


7.9- 10-* (1.0) 


0.011 (0.5) 



Table 6.3 
e = 0(10"^®), nonsmooth solution 



N 


\\u-Uh\\L^I) 


\W-Uh\\L^I) 


llk^llli 


IllellU 


100 


0.035 


0.070 


6.3-10-^ 


0.037 


200 


0.018 (1.0) 


0.048 (0.5) 


3.2-10-^ (1.0) 


0.024 (0.6) 


400 


9.1 • 10--^ (1.0) 


0.033 (0.5) 


1.6 • 10"^ (1.0) 


0.016 (0.6) 


800 


4.6-10-'' (1.0) 


0.023 (0.5) 


7.9- 10-* (1.0) 


0.010 (0.7) 



Table 6.4 
e = h, nonsmooth solution 
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